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The hydrogen plasma is studied in the very high density (atomic and metallic) regime by extensive 
ab initio Molecular Dynamics simulations. Protons are treated classically, and electrons in the 
Born- Oppenheimer framework, within the local density approximation (LDA) to density functional 
theory. Densities and temperatures studied fall within the strong coupling regime of the protons. 
We address the question of the validity of linear screening, and we find it to yield a reasonably 
good description up to r a w 0.5, but already too crude for r s = 1 (with r a = (S/Anp) 1 ^ 3 the ion 
sphere radius). These values are typical of Jovian planets interiors. Finite-size and Brillouin zone 
sampling effects in metallic systems are studied and shown to be very delicate also in the fluid 
(liquid metal) phase. We analyse the low-temperature phase diagram and the melting transition. 
A remarkably fast decrease of the melting temperature with decreasing density is found, up to a 
point when it becomes comparable to the Fermi temperature of the protons. The possible vicinity 
of a triple point bcc-hcp(fcc)-liquid is discussed in the region of r s ~ 1.1 and T ~ 100 — 200 K. 
The fluid phase is studied in detail for several temperatures. The structure of the fluid is found 
to be reminiscent of the underlying bcc (solid) phase. Proton-electron correlations show a weak 
temperature dependence, and proton-proton correlations exhibit a well-defined first coordination 
shell, thus characterizing fluid H in this regime as an atomic liquid. Diffusion coefficients are 
computed and compared to the values for the one-component plasma (OCP). Vibrational densities 
of states (VDOS) show a plasmon renormalization due to electron screening, and the presence of a 
plasmon-coupled single-particle mode up to very high temperatures. Collective modes are studied 
through dynamical structure factors. In close relationship with the VDOS, the simulations reveal the 
remarkable persistence of a weakly damped high-frequency ion acoustic mode, even under conditions 
of strong electron screening. The possibility of using this observation as a diagnostic for the plasma 
phase transition to the fluid molecular phase at lower densities is discussed. 
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I. INTRODUCTION AND STATE OF THE ART 

The possibility of hydrogen metallization under high 
pressure was first discussed by Wigner and Huntington 
in the thirties This particular subject was included 
in the more general conjecture that any system becomes 
metallic if sufficient pressure is applied. Electrons, which 
are bound in the isolated atomic or molecular species, 
hybridize in a condensed phase to form extended states 
which gather in energy bands. When pressure is ap- 
plied, the bands widen due to the enhancement of the 
overlap between neighboring atoms' orbitals, up to the 
point at which the energy gap between the valence and 
the conduction band vanishes, thus giving rise to metal- 
lic behavior. This enhancement of the overlap can also 
be viewed as the growing importance of the electronic 
kinetic energy relative to the ion-electron potential; the 
latter tends to bind electrons to the ions (also in the form 
of interatomic, intramolecular or intermolecular bonds). 
In the case of pure hydrogen, metallization was estimated 
to occur around 2.5 Mbars (1 Mbar = 100 GPa), i.e. at 
pressures that have become accessible to experiment in 
diamond anvil cells only very recently Q. The absence 
of clear signs of metallic behaviour up to 2.9 Mbar ||]4|] 
certainly adds to the challenge, and stimulates the con- 
tinuous improvement of experimental setups in the search 
for the reluctant metallization. 

An additional difficulty in the description of this par- 
ticular metal-insulator phase transition arises from the 
fact that at low pressures the low-temperature phase is 
not a monoatomic but a molecular solid, i.e. an hexago- 
nal close packed (hep) arrangement of H2 molecules M . 
In this sense, H behaves more like a halogen than like 
an alkali metal ||. Up to pressures of the order of 1 
Mbar, experiment indicates that the hep arrangement is 
preserved. At pressures of the order of 110 GPa, the 
low-pressure free-rotator phase of para-Eb freezes into 
an orientationally ordered phase, whose nature is not 
yet fully understood. The observation of more than 
one vibron mode Q precludes the hep structure with 
all the molecules pointing along the c-axis. However, 
other more complex hep structures with tilted molecules 
(herringbone- like) are compatible with experiment Q . It 
is also possible that orientational order is not complete, 
thus giving rise to a sequence of weak phase transitions 
until the perfectly ordered phase is reached. 

An intriguing phase transition at 150 GPa, signalled 
by a discontinuity in the molecular vibron Q|, has been 
ascribed either to a relative reorientation of the H2 
groups to metallization arising from an electronic 
band-overlap mechanism without molecular reorienta- 
tion H, and to both occuring at the same time Jl0| . 
Density functional (DFT) calculations indicate that ar- 
rangements with complex molecular orientations (non 
parallel) are energetically preferred in the high pressure 
phase ||. Another study |ll|] concentrated on other 
classes of structures including cubic Pa3 and rutile, and 



also two hep structures, namely Pca1\ and P1\jc. The 
hexagonal Pca2\ was found to be the favorable one above 
the 150 GPa transition. Very recently, Tse and Klug per- 
formed ab initio simulated annealing calculations with 96 
H atoms in the supercell, and found as a ground state 
an orthorombic structure, instead of hep, composed of 
groups of three strongly interacting if 2 molecules |L2j . 
All these candidates were found to be insulating, and 
the reason for this was claimed to be the opening of 
a hybridization gap at the Fermi level [[131 . Diffusion 
Monte Carlo (DMC) calculations by Natoli et al. El 
have essentially confirmed results by Kaxiras et al. W, 
also at the quantitative level for the angles between dif- 
ferent molecular units. They considered neither Nagara's 
Pca2i nor Tse's orthorombic structure. This DMC cal- 
culation also finds an insulating behaviour and, interest- 
ingly, zero-point-motion effects due to the protons appear 
to be quite structure independent. As a consequence, 
zero-point-motion would turn out to be irrelevant as re- 
gards the 150 GPa transition, contrarily to the claims of 
Surh et al. ffcj]. 

Well-converged local density functional (LDA) calcu- 
lations supplemented with the zero-point-motion energy 
contribution taken from frozen-phonon calculations [ p"5"[ , 
suggest that the molecular-hep phase goes over directly 
to an atomic squeezed-hexagonal phase, by breaking the 
intramolecular bonds, at a pressure of 380 GPa. How- 
ever, this calculation did not take into account non- 
hexagonal structures. Hence, also from the theoretical 
point of view, it is not clear up to now whether metalliza- 
tion happens already in the molecular phase, or whether 
it is a consequence of dissociation, thus occuring simul- 
taneously with the molecular-atomic transition A 
rhombohcdral structure was also postulated to supersede 
the squeezed-hexagonal phase at higher pressures (above 
400 GPa) p6[ , while the transition to a body-centered- 
cubic (bee) phase was located at around 1100 GPa [p^| . 

In general, static methods suffer from the drawback of 
having to carry out the investigation by guessing differ- 
ent structures. In this respect, the molecular dynamics 
results by Tse and Klug seem to be the most reliable. 
However, there is still some doubt because it is not obvi- 
ous that their 96-atom supercell with T-point sampling 
and fixed cell volume and shape calculation is a suffi- 
ciently good description. A recently proposed method 
that combines state-of-the-art electronic structure cal- 
culations with variable cell shape Molecular Dynamics 
simulations |l7f| , thus avoiding every undesired bias, is 
currently being used to investigate the low-temperature 
phase of H as a funcion of pressure jl8) . 

It is important to remark that, in general, energy 
differences between different structures are quite small, 
so that different levels of approximation often lead to 
different ground-state structures. In this respect, fully 
quantum-mechanical DMC calculations are very indica- 
tive. Ceperley and Alder jl9| located the molecular-to- 
atomic phase transiton at a pressure of about 300 GPa, 
by studying only cubic structures. Further refinements 
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of these DMC calculations have shown how dramatic the 
effect of improving the description can be pfj|] . Eventu- 
ally, it appears that the atomic ground state structure in 
the vicinity of the molecular-atomic transition is the di- 
amond structure, but the energy difference with respect 
to others is very small and, in particular, with respect to 
/3-Sn and hexagonal diamond, it is well within the error 
bars. The diamond structure was considered as a can- 
didate for the ground state in one of the previous DFT 
calculations, but found to be less favorable than a dis- 
torted hexagonal phase |D| . 

The monoatomic bec structure is unequivocally iden- 
tified as the ground state in the very high density limit. 
There, the electronic kinetic energy is largely dominant 
over the remaining contributions, such that the electrons 
behave essentially as free fermions; i.e. the electronic sub- 
system decouples from the protons and becomes a rigid, 
homogeneous electron gas, which acts only as a neutral- 
izing background for the protons. This is well-known as 
the one-component plasma (OCP) model, whose classical 
version was shown to crystallize into the bec structure. 
This follows from a simple calculation of the Madelung 
energy, although the energy differences relative to other 
structures turn out to be rather small. In the classical 
OCP, which was thoroughly studied during the past two 
decades the bec structure turns out to be the sta- 
ble one also at finite temperature and up to the melting 
point. The quantum version of the OCP at finite tem- 
perature has not yet been studied in detail. A numeri- 
cal study using the path integral Monte Carlo technique 
(PIMC) is currently under development |2^ |, 

At higher temperatures the protons begin to behave 
as classical particles. The degeneracy temperature for 
the protons can be estimated to be of the order of « 
180 K/r 2 , by comparing the mean interprotonic distance 
and the thermal de Broglie wavelength. The electron 
degeneracy temperature (T^) lies well above (by a fac- 
tor m p — 1836) such that, at temperatures lower than 
TJ w 326000 K/r 2 s — 1 Hartree/rg, electrons can be con- 
sidered to follow the protons adiabatically, always staying 
in the the ground state compatible with the current pro- 
tonic configuration (Fermi temperatures are about twice 
these numbers, e.g. T F = 326 K/r 2 ). This latter ap- 
proximation allowed Hohl et al. [E3| to examine the hot, 
dense molecular phase by means of ab initio Molecular 
Dynamics (AIMD) simulations, and also to make some 
progress regarding the low-temperature structure. The 
AIMD study of the atomic phase is the subject of this 
and of a previous publication . The effect of excited 
electronic states in the same regime has very recently 
been explored by a novel ab initio Molecular Dynamics 
scheme using Mermin's density functional [ p5| , and the 
fully-quantum PIMC method was applied to simulate the 
very high-temperature fluid phase |26| , where the number 
of excited electronic states involved becomes too large to 
be treated efficiently with Mermin's functional. The low- 
temperature regime, where protons in turn become de- 
generate, requires different simulation techniques which 



are currently being developed |27[ |. 

High-temperature high-pressure measurements on Hy- 
drogen have been very recently reported in shock wave 
experiments pq| . Pressures between I and 2 Mbar and 
temperatures of some thousand K are now feasible and 
rather controllable, and further expansion of this range 
is in sight. In such experiments, Nellis et al. observed 
the metallization of fluid molecular Hydrogen at P = 1.4 
Mbar and T = 3000 K, a temperature which is signif- 
icantly lower than that predicted for the plasma phase 
transition by Chabrier et al. j29| on the basis of a multi- 
component thermodynamic theory that employs an ap- 
proximate equation of state. Although more accurate 
PIMC simulations have essentially confirmed the predic- 
tions of the theory, including the existence of a first order 
phase transition at high pressures and a somewhat lower 
value for the transition temeprature [fj0| , the discrepancy 
with respect to shock wave experiments and AIMD re- 
sults |31j still holds. The latter would rather indicate a 
continuous transition in which metallization and disso- 
ciation of H2 molecules are closely related phenomena. 
A possible scenario would be the existence of two phase 
transitions: a continuous one, at relatively low tempera- 
tures, where molecular Hydrogen already metallizes and 
begins to dissociate, and a discontinuous one at higher 
temperatures where a massive dissociation occurs with a 
concomitant jump in the electrical conductivity. 

Hydrogen at finite temperatures and high densites con- 
situtes, hence, a strongly coupled proton-electron plasma 
which is of great astrophysical interest since it is, in par- 
ticular, the major constituent in Jovian planet interiors. 
Actually, the latter are basically H plasmas with a few 
percent admixture of He, and the statistical properties 
of the mixture (mixing/demixing transition) are of great 
importance to account for experimental observations and 
to establish models for the evolution of these planets |}2| . 
Nevertheless, the answer to many relevant questions in 
astrophysics requires a full understanding of the statisti- 
cal properties of the pure main constituents. Hydrogen is 
present in both atomic and molecular fluid phases. These 
are separated by a boundary located at some distance 
from the center of the planet. The characteristics of this 
boundary, e.g. precise location, width, etc. depend on 
the density and temperature profile, i.e. on the equation 
of state. The matter of dissociation and metallization 
is particularly relevant because the large magnetic field 
measured in Jupiter would include a non neglgible com- 
ponent generated in the outer molecular fluid phase, pro- 
vided that temperature is above the metallization thresh- 
old. To quantify this effect the important quantities are 
the dissociated fraction and the electrical conductivity, 
which can be obtained from simulations in the molecular 
phase. 

The hydrogen plasma is perhaps the most fundamen- 
tal, simple many-body system, with all the interactions 
(proton-proton, proton-electron and electron-electron) 
described exactly by the bare Coulomb potential. In spite 
of this, the phase diagram appears to be surprisingly rich. 
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The purpose of this article is to present a detailed study 
of the very high pressure physical properties and phase 
diagram of this simple and fascinating system. We are 
going to deal always with atomic phases, the molecular 
phases appearing at lower densities than those studied 
here. In Section II we briefly describe some details of 
the simulations, and we address some important techni- 
cal issues such as the interplay between pseudopotentials 
and basis set expansions. Section III is devoted to the 
analysis of the validity of linear response theory in its 
ability to describe proton-proton interactions in terms of 
an effective pair potential of the screened Coulomb form. 
In Section IV we address a crucial problem that arises in 
the ab initio simulation of metallic systems, in particular 
liquid metals, i.e. size effects and Fermi surface sampling. 
In Section V we enter directly into the properties of the 
Hydrogen plasma by describing the solid phases and the 
melting at very high pressure. Section VI is concerned 
with the fluid phase, which we characterize as atomic- 
like. Diffusion and vibrational properties are presented 
in Section VII, while Section VIII is devoted to a thor- 
ough study of the collective dynamics, as a approach to 
the metal-insulator transition coming from the metallic 
side. The behavior of the longitudinal collective mode 
is studied as a function of density and temperature, and 
proposed as a probe for the metal-insulator transition at 
finite temeprature. Finally, we conclude in Section IX. In 
the remaining part of this Section we define some useful 
parameters and ratios. 

The plasma is made up of n electrons and as many 
ions (protons) per unit volume, such that the usual 
dimensionless density parameter is r s = a/as = 
(3/47rn) 1 / 3 /as, where a is the mean interionic distance 
(ion-sphere radius) and as is the Bohr radius. Adopt- 
ing atomic units throughout, the Fermi momentum is 
kp = (9/47r) 1 / 3 /r s , and the Thomas-Fermi screening 
length is kpF = (12/n) 1 / 3 / y/r^. The dimensionless 
Coulomb coupling constant associated with the classical 
ions is T = e 2 /(ak B T). Note that k B T = l/(Tr s ), and 
that the electron degeneracy parameter is 9 = T/Tp — 
2(4/97r) 2 / 3 r s /r. Typical densities inside Jovian planets 
are between 1 and 10 g/cm 3 and temperatures are of 
the order of 100 to 10000 K, i.e. r s ranging from 0.5 to 
1.5, and T w 20 - 200 Note that within this range 
9 w 0.015 <C 1 and for the present calculations, whose 
aim is precisely to address some aspects of astrophysi- 
cal plasmas, we can resort to the adiabatic approxima- 
tion by assuming that the electrons are always in their 
instantaneous ground state for any given ionic configura- 
tion. Higher temperatures require the relaxation of this 
hypothesis. 

II. TECHNICAL DETAILS OF THE 
SIMULATIONS 

Our simulations were carried out using the stan- 
dard Car-Parrinello (CP) AIMD scheme f|. Thc 



ions were considered as classical point-like particles, 
whereas the electrons were described by density func- 
tional theory (DFT) within the local density approxima- 
tion (LDA) (34j. The electronic density was constructed 
via Kohn-Sham single-electron orbitals, expanded in a 
plane wave basis set up to a specified energy cutoff -E C ut 
(see below). We have extensively studied a system con- 
sisting of 54 H atoms contained in a simple cubic sim- 
ulation box under periodic boundary conditions (PBC), 
and we have also performed simulations for larger su- 
percells, containing 128 and 162 H atoms, in order to 
analyse finite-size effects. The plane wave expansion was 
carried out around the T-point of the supercell's Brillouin 
zone. Finite-size and Brillouin zone sampling effects will 
be analysed later, in a forthcoming section. 

We used the exchange-correlation functional deduced 
from the results of quantum Monte Carlo calculations on 
the uniform electron gas |56| , and the bare Coulomb po- 
tential for the ion-electron interaction. The singularity 
of the ion-electron Coulomb attraction implies that there 
is never absolute convergence of the electronic quanti- 
ties as the energy cutoff for plane waves is increased. In 
particular, the cusp condition is never satisfied, the den- 
sity reaching the origin with zero slope. This is because 
the Fourier expansion of 1/r is also a long-range func- 
tion (47r/(? 2 ), meaning that there always exists a spa- 
tial scale small enough to require a representation in- 
volving large-<? components. To satisfy the cusp condi- 
tion, a different (localized) basis set is needed instead of 
plane waves, e.g. Slater-type orbitals |37|| . A finite cutoff 
g cu t = ^2E cut (Ry) translates into the following effective 
potential: 

„„„„.) _!(,_> r iiw dr \ (1) 

r \ * J gcut x J 

This description is variational in the number of plane 
waves (i.e. in the energy cutoff) and, even if the energy 
never fully converges, other quantities, like the electronic 
distribution around the protons, are affected only in a 
tiny region around the nuclei, the remainder being es- 
sentially converged at a finite cutoff. This is shown in 
Fig. 1 for protons fixed at the ideal bec lattice sites. The 
choice of the cutoff strongly depends on the average den- 
sity: the lower r Sl the higher the cutoff. It is the number 
of plane waves which has to be kept approximately con- 
stant in order to achieve similar levels of convergence. In 
order to have the electronic density correctly described 
from less than one third of the nearest-neighbor ionic 
distance onward, we decided to vary the cutoff between 
230 Ry for r s = 0.5 and 60 Ry for r s = 1 and 1.2. It 
can be seen in Fig. 1 (r s = 0.5) that differences in the 
proton-electron correlation function become noticeable at 
distances smaller than 0.5 a, while the nearest neighbor 
distance in the bec lattice is d nn = 1.76 a. It should be 
pointed out that a cutoff of 60 Ry is sufficient to achieve 
convergence in the properties of the H2 molecule, like 
the proton-proton distance (1.5 a.u. = 0.79 A) and the 
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vibrational frequency (4160 cm" 1 ). The difference with 
respect to the experimental quantities (0.74 A and 4400 
cm -1 ) may be ascribed exclusively to the LDA. 

Isothermal ionic dynamics (particularly in the fluid 
phase) was achieved by using a Nose-Hoover thermostat. 
The time step for integration of the CP equations of mo- 
tion was chosen typically between 0.25 and 1.5 atomic 
units (depending on the temperature and density), i.e. 
roughly 10 -17 s, compared to the ionic plasma period 

3/2 

Tp « 155 ry a.u. In some cases another Nose thermo- 
stat was necessary to keep the fictitious kinetic energy 
of the electronic degrees of freedom at low values. Most 
runs extended over 0.4 to 0.6 psec, which amounts to 
more than 50 plasma oscillations. Before this, we al- 
lowed for thermal equilibration during an initial period 
of 0.2 psec in the solid phase, and 0.8 to 1.2 psec in 
the fluid phase. The densities r s — 0.5, r s = 1 and 
r s = 1.2 were explored in detail for temperatures in the 
range 100 K < T < 10000 K, which correspond to the 
strong coupling regime (T > 30), and cover both fluid 
and solid phases of the ionic component. In particular, 
r s 1 and T w 7000 K are typical conditions inside 
Jovian planets |3S{| . The simulations at r s — 1.2 were 
carried out only in the fluid phase because the stable 
solid has no longer the bec symmetry, and thus it is not 
compatible with the choice of 54 H atoms in a simple 
cubic supercell. 



III. VALIDITY OF LINEAR RESPONSE THEORY 

In the very high density limit (r s — > 0) the electrons 
are barely polarized by the ionic charge distribution, due 
to their rapidly increasing Fermi energy, and the hydro- 
gen plasma practically reduces to two decoupled compo- 
nents: a classical ionic "one-component-plasma" (OCP), 
and a degenerate, rigid electron "jellium" . Both systems 
have been extensively studied by classical and quan- 
tum |3q| computer simulations. For finite, but small r s 
(r s <C 1), the coupling between the two components 
may be treated perturbatively within linear screening 
theory [^9|. However, due to the strength of the un- 
screened Coulomb interaction between protons and elec- 
trons, such a perturbative approach is expected to break 
down rapidly as r s increases. In this section we investi- 
gate the high density limit and the validity of the descrip- 
tion of the hydrogen plasma in terms of a linear response 
picture (LRT) for the electronic component. 

To this end we compare the proton-electron pair corre- 
lation function taking into account the full LDA response: 
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and the one obtained by replacing the electronic charge 
density by its linear response expression in terms of the 
ionic charge density: 



p e (k)= X e(k)[- 
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(3) 



where Xe( k ) is the static electron density response func- 
tion. In eq. (2), jo(x) denotes the I = spherical Bessel 
function, pi (k) is a Fourier component of the microscopic 
ionic density, and p e (— k) = /9*(k) is a Fourier component 
of the electronic density corresponding to the instanta- 
neous ionic configuration {R/}. The latter represents, 
within the adiabatic approximation, a ground-state ex- 
pectation value. 

For the density response we adopt the random phase 
approximation (RPA) corrected by the long wavelength 
limit of the local field function G(k) pl| . 
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with l(x) the Lindhard function and kp F the Thomas- 
Fermi wavenumber. For the local field correction we 
use the following expression, which is consistent with the 
LDA and Slater's local exchange: 



G(k/k F ) 
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where the correlation energy density e c {r s ) is that pro- 
posed by Perdew and Zunger Q, and A = (A/9tt) 1/3 . At 
typical densities studied in this work the exchange con- 
tribution to G{k/kp) turns out to be largely dominant 
over correlation, and both of them represent a small cor- 
rection relative to the RPA susceptibility. In fact, this 
is reasonable because the RPA (or Lindhard approxima- 
tion) is known to be good in the high density limit. 

The proton-electron distribution function reads ac- 
cordingly: 
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In Fig. 2 we compare the proton-electron radial dis- 
tribution functions in the ideal bec structure as obtained 
from the LDA and from LRT, for r s — 0.5 and r s = 1. 
As expected, g pe (r) is considerably more structured at 
the lower density. The LDA and LRT distribution func- 
tions are very close at r s = 0.5, but significant differences 
are clearly apparent at r s — 1, reaching a region well be- 
yond the first ionic shell and signalling the break-down of 
the linear screening regime. When the ions are at finite 
temperature the differences are significantly enhanced, 
particularly at short distances. Also the location of the 
first minimum turns out to be shifted outwards by LRT 
in the fluid phase (unlike in the solid phase), thus indi- 
cating that the LRT description of the electronic charge 
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distribution worsens for increasing temperature. This is 
shown in Fig. 3. The dashed curves have been obtained 
by averaging the adiabatic electronic charge distribution 
around the protons along the AIMD trajectory, assum- 
ing that the configurations generated in this way are also 
representative of the LRT. In fact, the magnitude of the 
differences observed in Fig. 3 suggests that the LRT 
trajectories will differ significantly from the LDA ones, 
implying that this averaged electronic distribution might 
be meaningless within LRT. In other words, the fully- 
consistent LRT is likely to be worse than the rasults pre- 
sented here. 

Following a novel procedure |0J we have fitted, for 
r s = 1, an effective proton- proton pair potential to our 
set of AIMD configurations generated at several temper- 
atures spanning both solid and fluid phases. In Fig. 4 we 
compare this potential to the one obtained within LRT, 
i.e. by Fourier transforming ulrt (k) = 4ir[l — Xe(k)]/k 2 , 
where Xe(k) is given by expression (Q). The nearest 
neighbor distance in the bcc structure for r s = 1 is 1.76 
a.u., and the position of the first peak in the proton- 
proton radial distribution function decreases down to val- 
ues of the order of 1.5 a.u. in the fluid phase and upon 
heating (see below). This means that nearest pairs of 
protons will spend most of the time at distances of this 
order. Looking carefully at Fig. 4 it can be seen that, at 
those distances, the two potentials differ by more than 
10%, the LDA one being steeper. On the other side, the 
LRT pair potential appears to be more long-ranged than 
the LDA one. This implies that the LDA pair potential 
has a characteristic screening length shorter than its LRT 
counterpart, i.e. screening is more efficient than linear at 
r s = 1. 

The departure from the linear screening regime, also 
in the form of many-body effects beyond the pair po- 
tential approach (e.g. three-body terms or embedding 
functions), becomes more pronounced as r s is increased, 
leading to ground state-atomic phases other than bcc, fee 
or hep (e.g. hexagonal and diamond) and, eventually, to 
recombination and to the different H2 molecular phases. 
The reason for this early departure from the LR regime 
has to be traced back to the unusual strength of the bare 
Coulomb interaction, arising from the absence of core 
electrons. In fact, other alkali metals can be reasonably 
well described by LRT at much lower electronic densities 
(r a > 3) 42j. Interestingly, the range of validity deduced 
here for the LRT is much wider than expected in previous 
theoretical work based on perturbative expansions [j39| , 
where r s — 0.1 was identified as the upper limit. 



IV. SIZE EFFECTS AND FERMI SURFACE 
SAMPLING 

In the solid phase, the combination of T-point sampling 
and a finite system size is not expected to provide a very 
accurate description of the electronic component, because 



all the periodicities beyond the size of the supercell are 
not properly included. This is particularly important in 
metallic systems, where no point in the Brillouin zone can 
be taken as representative of the band structure of the 
bulk solid. The reason is that occupied states (contribut- 
ing to the electronic density) and empty states (which 
do not contribute) coexist in the same band, correspond- 
ing to different k-points. A particular choice (e.g. the 
T-point and its refolded images) will sample the conduc- 
tion band in some specific points, but the character of the 
Fermi surface could be misinterpreted if empty states are 
taken as if they were occupied and viceversa. For quasi- 
spherical Fermi surfaces this is unlikely to happen, but 
for transition metals or semimetals (like graphite) this 
effect is crucial, and a very fine sampling of the Brillouin 
zone is needed to obtain the right physics. 

In the case of simple metals the effect of quantization of 
the electronic states in the simulation box is more impor- 
tant. The Brillouin zone of the unit cell is sampled with 
a finite number of points, which arise from the refold- 
ing of the T-point of the supercell. These points reflect 
the symmetries of the system in the sense that, if the 
T-point refolds onto some point ko, then all the points in 
the star of ko must arise from some other refolding of T. 
Otherwise, the symmetry is broken and spurious forces 
appear that drive the system away from the symmetric 
ground state. Depending on the case at hand, the dis- 
tortion can be rather large, especially if the symmetry of 
the supercell is very different from that of the unit cell. If 
the correct symmetry is used for the supercell then, since 
all the points in the star are equivalent, the eigenvalues 
associated with them are degenerate, thus giving rise to 
the formation of electronic shells. If a better sampling 
of the Brillouin zone is performed, or the size of the sys- 
tem is increased, new shells appear that eventually give 
rise to a continuous energy band. But for finite systems 
and restriction to the r-point, the electronic density of 
states consists of a set of discrete peaks (the shells). A 
problem arises with the highest occupied shell because 
in finite-size metallic systems it is partially occupied, un- 
less a very fortuitous situation occurs. The immediate 
consequence is that it is not possible to fullfil the sym- 
metry requirements by occupying all the points in the 
star with an integer number of electrons. Then, unless 
fractional occupation is introduced, the symmetry is bro- 
ken even when the symmetry of the unit cell is mantained 
for the supercell. This is analogous to the Jahn- Teller ef- 
fect in molecular systems, where a lower-energy state can 
be obtained by reducing the symmetry and breaking the 
electronic degeneracy. However, in the present case this 
effect is unphysical. 

In a fluid phase these effects are normally expected to 
become less important because of the breakdown of the 
discrete translational invariance, as reflected in the very 
existence of a finite Brillouin zone via Bloch's theorem. 
The unit cell becomes of infinite size and the Brillouin 
zone reduces to a single point, i.e. the F-point. However, 
for practical reasons, computer simulations are bound to 
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represent the infinite (fluid) system with a limited num- 
ber of particles (typically of the order of some hundreds 
in ab initio calculations) , while repeating periodically the 
supercell (via PBC). In fact, this introduces a spurious 
Brillouin zone, which is associated with periodicities that 
are absent in the infinite fluid. The properties of the fluid 
are assumed to be recovered from the finite sample in 
terms of PBC combined with statistical averages, which 
in the case of our AIMD simulations are computed as 
time averages. 

This is a reasonable justification for a purely classical 
system and also for liquid semiconductors but, for metal- 
lic systems, as soon as electronic states are introduced, 
the quantization of these states in the (unphysical) simu- 
lation box acquires a crucial role. The problem of partial 
occupation of the star persists in the fluid, in the sense 
that spurious forces appear that modify in a non-trivial 
way the structural and dynamical properties of the sys- 
tem. 

A qualitative picture of the consequences of these ob- 
servations can be obtained by comparing the proton- 
proton radial distribution function for systems of differ- 
ent sizes. A system of 54 H atoms (of valence charge 
equal to 1) is really a fortuitous case of compatibility of 
a closed electronic shell (the whole star is fully occupied) 
with a bcc atomic arrangement in a simple cubic super- 
cell. In the case of 128 H atoms only 57 states are doubly 
occupied, and the next shell of 24 degenerate states has 
to accomodate 14 electrons (i.e. only 30% of the shell). 
On the contrary, 162 H atoms is not compatible with a 
bcc arrangement in a simple cubic supercell, but it closes 
the former electronic shell thus amounting to 81 doubly 
occupied states. Therefore, in the fluid phase it should 
behave essentially as the 54-atom supercell but better 
converged in system size (or k-points). In Fig. 5(a) we 
show the curves corresponding to these 3 different sys- 
tem sizes at r a = 1 and T — 1000 K. In fact, up to the 
boundary of the 54-atom cell (r — 3 a.u.) there is hardly 
any relevant difference between the g PP (r) correspond- 
ing to 54 and 162-atom supercells. The same kind of 
picture holds also at higher temperatures, meaning that, 
as regards static thermodynamic properties, 54 H atoms 
already give a very reasonable picture. 

Very different is the situation with open shell systems 
like the 128 H atoms one (short-dashed line in Fig. 5(a)). 
It is clear that this g PP (r) has little to do with those of 
54 and 162 atoms. The first peak is significantly lower 
and broader, and also its position is shifted downwards. 
The first valley is much shallower and, in practice, the 
pair distribution becomes structureless beyond it, while 
with 162 atoms the atomic shell structure is still visible at 
least up to the second valley. The 128-atom distribution 
resembles, in fact, the one that would have been obtained 
with closed shells at a higher temperature value. 

The correct way to account for open shell structures 
at finite temperature is to doubly occupy the lowest- 
lying electronic states at every time step of the AIMD 
simulation, while keeping empty the rest of that shell. 



The curves in Fig. 5(a) have been obtained by apply- 
ing the standard procedure of considering explicitly, in 
the description of the electronic component, a number of 
electronic states equal to one half of the number of elec- 
trons, i.e. including those strictly necessary. The justifi- 
cation for this approach is that typical temperatures are 
much lower than the Fermi temperature. Atomic mo- 
tion leads to a degeneracy lifting of the order of a few 
percent of an eV, implying that the distance between 
electronic states is usually much larger than the width of 
the Fermi-Dirac distribution; this latter is of the order 
of k B T, i.e. about 0.1 eV at T = 1000 K. This leads 
to a symmetry breaking in the sampling of the Fermi 
surface, which is expected to be recovered in terms of 
statistical averaging. However, following the approach 
above, in practical simulation times we have not noticed 
a convergence to the closed-shell picture. In Fig. 6 we 
show the Fermi-Dirac distribution function for three dif- 
ferent temperatures, namely 1000 K, 5000 K, and 10000 
K, compared to the electronic density of states averaged 
along an AIMD run at T — 1000 K, and to a snapshot 
of the instantaneous Kohn-Sham eigenvalues. It is clear 
that the occupation numbers fall from 2 to in an energy 
scale narrower than the thermal splitting of the eigenval- 
ues. 

The drawback of this kind of approach is that for 
open shells the ordering of the states within the high- 
est (partially) occupied shell changes continuously during 
the MD evolution and, in particular, occupied states be- 
come empty and viceversa. The standard CP approach 
is not able to take into account this phenomenon, and 
this is one reason for the well-known failure in metallic 
systems, reflected in the energy transfer between elec- 
tronic and ionic degrees of freedom (cooling the ions and 
heating the electronic orbitals). The double Nose ther- 
mostat proposed by Blochl and Parrinello |43| helps in 
fixing the temperature of each of the components, but 
still does not take into account level crossing effects. To 
our knowledge, there are three methods capable of solv- 
ing this problem: one is to abandon the CP Lagrangian 
strategy in favor of a self-consistent minimization at each 
time step Q , a second one is to abandon the description 
of the electronic component in terms of single-particle or- 
bitals to work directly with the electronic density which 
is, by definition, constructed with the lowest occupied 
states [ p5| , and finally a third one consisting of a rigor- 
ous Lagrangian formulation which incorporates the oc- 
cupation numbers (actually the conjugate variables, i.e. 
the Kohn-Sham eigenvalues) as dynamical variables Es[ . 
However, this latter procedure exhibits the odd feature 
that the fictitious kinetic energy of the electronic orbitals 
still increases at the expense of the ionic component, due 
to the appearance of low-energy excitations introduced 
by the dynamics of the eigenvalues. 

An approximate solution can still be found within the 
Lagrangian approach by evolving explicitly the whole 
open shell, but initially occupying - with an integer num- 
ber of electrons - only the lowest states of the this shell. 
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This will approximately take into account level crossing 
in terms of mixing of states within the shell, during the 
time evolution. The very same mechanism that leads to 
energy transfer between ions and electronic orbitals, i.e. 
the vanishing energy gap Eg ], is responsible for mixing 
occupied and empty orbitals, thus allowing for initially 
empty states (if explicitly included) to become occupied 
and viceversa. Statistical averaging completes the task 
by generating a more uniform sampling of the Fermi sur- 
face. Fig. 5(b) shows how the results obtained with 
128 H atoms can reproduce approximately those obtained 
with 54 and 162 H atoms. It has to be pointed out that 
this procedure, although not rigorously justified, has the 
nice feature that the fictitious kinetic energy of the elec- 
tronic orbitals is practically constant, behaving exactly 
as in systems with a gap, i.e. there is no energy exchange 
between ionic and electronic degrees of freedom. 



V. LOW-TEMPERATURE PHASE: STRUCTURE 
AND MELTING 

In the OCP (r s — > 0) limit the ionic bcc lattice is 
known to be the stable crystalline structure up to melt- 
ing (which occurs at T w 180 |21|). We have studied the 



stability of the bcc structure at low temperatures and fi- 
nite r s by performing canonical MD simulations; initial 
conditions were constructed by giving the ions a small 
random displacement (about 3 % of the nearest neigh- 
bor distance d nn ) from their alleged equilibrium positions 
(the bcc lattice sites). The bcc structure was found to 
be dynamically stable against such displacements at least 
up to r s = 0.5. At r s = 1 the bcc structure was found to 
be unstable for T < 100 K, where a close-packed struc- 
ture appears to be favored. This is consistent with a de- 
scription in terms of LRT. In fact, at high densities the 
effective LRT potential behaves essentially like a Yukawa 
potential 



vsc(r) = -exp(-r fc TF ) 
r 



(7) 



while Friedel oscillations are practically negligible. The 
phase diagram of a classical system of particles inter- 
acting via the above Yukawa potential (Q) has been ex- 
tensively studied by computer simulations and lattice 
dynamics [^tJ. These calculations point to a bcc-fcc 
phase transition when the density-dependent screening 
wavenumber increases, i.e. when the effective interaction 
becomes of shorter range. At T = the bcc structure 
is found to be stable up to r s 0.6, beyond which the 
fee phase becomes the stable structure. The r s at co- 
existence shifts to higher values at finite temperatures, 
such that the system goes through a structural fcc-bcc 
phase transition as the temperature increases along an 
isochore. 

The situation here is reminiscent of the behavior of 
alkali metals. Na exhibits an hep ground state, while Li 



goes from bcc to fee and eventually to hep at very low T. 
The heavier alkalis K, Rb and Cs undergo a structural 
transition to fee upon cooling below T w 5 K . In all 
these cases the entropic contribution of the bcc structure, 
arising from the valley in the phonon dispersion along 
the (110) direction, wins over at finite temperature and 
stabilizes this phase. This is exactly what we observe in 
our simulations for H, where at r a = 1 the bcc structure 
appears to be stable for T > 100 K. 

Still, finite energy cutoff, finite system size, and coarse 
Brillouin zone sampling may have a large influence on 
the stability of the ground state structure. The study 
of this part of the phase diagram deserves special atten- 
tion because also zero-point motion effects on the pro- 
tons have been shown to influence the stability of differ- 
ent structures at T = (2(J. Disregarding the problem 
of zero-point energy (ZPE), and only as a check of the 
present calculations (which do not include ZPE), we have 
performed total energy full-potential Linear Muffin Tin 
Orbitals (FP-LMTO) calculations for solid, monoatomic 
H in the bcc, fee and hep structures, at r s — 0.5 and 
r s = 1. The energy differences turned out to be very 
small, but at r s — 1 the fee and hep structures are sig- 
nificantly lower in energy than the bcc. Moreover, the 
energy differences are enhanced if a coarse sampling of 
the Brillouin zone is performed. These calculations also 
identify the hep structure as the lowest energy one, but 
the difference with respect to the fee is within the accu- 
racy of the calculations. The reason for this can be found 
in Fig. 4. The differences between fee and hep structures 
begin only at the level of third nearest neighbors, a re- 
gion where the effective potential shows negligibly small 
Friedel oscillations. 

Between r s = 1 and 1.2 a phase transition occurs that 
takes the system from hep to a simple-hexagonal phase 
with a compressed c/a ratio (squeezed-hexagonal). This 
is compatible with static total energy calculations by Bar- 
bee et al. ]l5j , and it is an additional confirmation of 
the breakdown of LRT, because pair potentials (like the 
LRT one) are not able to stabilize anisotropic structures 
like the simple-hexagonal. A more detailed study of the 
low-temperature atomic phases of Hydrogen is currently 
under way Q . 

Next, the melting of the ionic crystal was investigated 
by gradually increasing the temperature and monitor- 
ing the time-dependent mean square displacement of the 
ions <| Ar(i) | 2 >=<| r(i) - r(0) | 2 >. In the crystalline 
phase, <| Ar(t) | 2 > goes over to 2 <| Sr | 2 > for suffi- 
ciently long times, where <| <5r | 2 >=<| r — R | 2 > de- 
notes the static mean square displacement (R are the 
equilibrium positions of the ions). In the fluid phase dif- 
fusion sets in, so that <| Ar(i) | 2 >= 6Dt at long times, 
with D the ionic diffusion constant. At r s = 0.5 dif- 
fusion was found to set in at T 230, which may be 
identified with the limit of mechanical stability of the 
(overheated) metastable crystal. The thermodynamic 
transition occurs at lower temperature (higher T); its 
location may be estimated by assuming that the Lin- 
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demann ratio L — (<| Sr | 2 >) 1 / 2 /c? at melting is the 
same as for the OCP, i.e. L w 0.15 j3f|. This leads to 
T m (r s = 0.5) « 290 compared to T m (r s = Q) w 180, 
indicating a strong influence of electron screening on 
the melting transition. This confirms recent predictions 
based on free energy comparisons, obtained by means of 
an approximate density functional theory p5[ . 

We have also studied the melting transition at r s = 1 
using the same procedure; the melting temperature drops 
sharply from T m (r s = 0.5) w 2200 K to T TO (r s = 1) w 
350 if (r m (r s = 1) ps 930). This indicates the possible 
proximity of a triple point bcc-hcp(fcc)-liquid, analogous 
to that found for Yukawa potentials. However, the non- 
linearity of the screening at these values of r s is likely 
to bring the triple point from r s = 3 |47| to r s ps 1.1. 
Moreover, the hep structure goes over to a simple hexag- 
onal one at r s somewhere between 1 and 1.2, and this 
has to be a consequence of the appearance of anisotropic 
forces, beyond the level of pair-wise additivity. The very 
low value of the melting temperature might also be re- 
lated to the appearance of these forces. Interestingly, at 
r s = 1 the Fermi temperature of the ionic component is 
Tp = 326 K, a value close to the melting point (350 K). 
Therefore, the influence of quantum effects for the pro- 
tons on the melting transition cannot be ignored, and 
will probably also play an important role in the above 
structural phase transition. In particular, they might 
destabilize the hexagonal phase in favor of some more 
isotropic configuration which has a higher energy within 
a framework of classical protons. 

VI. THE FLUID PHASE: AN ATOMIC-LIKE 
PLASMA 

Turning to the fluid phase, a quantitative measure 
of ion-electron correlations is provided by the spherical- 
ized average of the ion-electron pair correlation function 
gpe(r), as computed from eq. (2). The effect of tem- 
perature on g pe (r) is illustrated in Fig 7, for r s = 0.5 
and r s = 1. The distribution function is seen to be re- 
markably insensitive to T over the whole range of tem- 
peratures, covering the solid and fluid phases, as already 
noticed at lower T (higher T) by Dharma-Wardana and 
Perrot in the framework of an approximate static DFT- 
HNC calculation |4§]. The observed weak temperature 
dependence implies that the main effect of ionic thermal 
motion (electrons are always at T = 0), is to enhance the 
localization of the electronic charge close to the protons. 

This behavior is to be contrasted with the predic- 
tions of higher-level theories that go beyond the Born- 
Oppenheimer approximation by including excited elec- 
tronic states. Both, fully-quantum PIMC Q| and Mer- 
min functional [ p0[ simulations imply that ion-electron 
correlations become weaker as temperature increases. 
This is to be intuitively expected, but for temperatures 
much higher than the ones studied here. Excited elec- 



tronic states are much more insensitive to ionic polar- 
ization effects, because they correspond to larger kinetic 
energies. In this way, for temperatures larger than a 
threshold value that can be estimated around 9 = 0.1, 
i.e. T s» 60000 K for r s = 1, the above localization ef- 
fect due to ionic disorder starts to be compensated by 
the effect of electronic excitations, so that eventually the 
opposite trend will take over. 

The ion-electron pair correlation functions for all tem- 
peratures are seen to intersect at a well-defined (reduced) 
distance from the proton site, irrespective of thermal 
ionic disorder and almost independently of density. We 
locate this value at r* « 1.3 a, and notice that the ratio 
t* /d nn w 0.73 is related to the ratio of the location of 
the nodes corresponding to the first two spherical Bessel 
functions (jo(x) and ji(x)). In fact, the electronic prob- 
lem can be modelled, in a very crude approximation, as 
that of a particle in a spherical well; the corresponding 
radial solutions are precisely the spherical Bessel func- 
tions. Temperature effects can be mimicked by increas- 
ing the relative population of excited states with respect 
to jo(x). However, the location of the node of ji(x) (the 
leading excitation) relative to jo{x) (the ground state) 
does not depend on temperature. Since the location of 
the nodes is defined in units of the radius of the well, and 
this is identified with d nn , which is proportional to r s , the 
crossing should not depend significantly on density. The 
first maximum of g pe (r) is clearly associated with the lo- 
cation of the first coordination shell (the first maximum 
in g pp (r) — see below), which is quite natural since the 
electronic density peaks at the proton sites. It is interest- 
ing to notice that differences in the electronic screening 
properties between r s = 1 and r s = 1.2 are significative 
only in the vicinity of the protons, up to r rs 0.5 a. 

The proton-proton pair correlation function for r g = 1 
is illustrated in Fig. 8(a), as a function of tempera- 
ture. It can be observed that the first peak remains 
clamped at the nearest neighbor distance (d nn = 1.76) 
for temperatures below the mechanical stability limit 
T s (r s = 1) rs 500A" -, while in the fluid phase it shifts 
continuously to shorter distances. The same plot shows 
that the location of the first minimum is quite insensitive 
to temperature. This, together with the fact that equiva- 
lent results are found at r s = 0.5, defines quite univocally 
a first coordination shell of radius rf cs w 2.4 a (with a 
the ion-sphere radius). The integrated number of parti- 
cles is shown in Fig. 8(b) as a function of temperature. 
The main result is that the first coordination shell con- 
tains 14 atoms on average, implying that the short-range 
structure of the liquid is quite reminiscent of that of the 
solid, since the first coordination shell of the bec struc- 
ture (containing first and second nearest neighbors) also 
contains 14 atoms. Simulations performed with 162 H 
atoms show that this is a genuine feature and not an 
artifact of the small system size. A second coordina- 
tion shell is also well-defined in the fluid phase provided 
that the temperature is low enough, i.e. T < 1500 K 
(r > 200), as can be observed in Fig. 9. However, the 







fluid becomes structureless beyond the first coordination 
shell at temperatures of the order of 5000 K (at r s = 1). 
Summarizing, the fluid phase of the H plasma at moder- 
ately high temperatures and very high densities (typical 
of the inner H shell of Jovian planets) behaves like an 
atomic liquid with a well-defined first coordination shell. 

The influence of electron screening is also apparent 
when comparing the ion-ion and charge-charge static 
structure factors Su(k) and Szzik). While at r s = 0.5, 
the two are nearly indistinguishable, the amplitudes of 
their main peaks differ significantly for r s = 1 (by roughly 
9 %), and r s = 1.2 (12%). Due to the discrete sampling of 
the electronic density in reciprocal space, the curves are 
noisy and will not be reproduced here. The strong polar- 
ization of the electronic component leads to a damping of 
the local charge fluctuations, and hence to a reduction of 
Szz(k). Again, the presence of a well-defined first peak 
and valley in S(k) is an indication that the fluid is well 
structured in this region of the parameter space. 



VII. DIFFUSION COEFFICIENTS AND 
VIBRATIONAL PROPERTIES 

Our AIMD simulations give direct access to the ionic 
dynamics by analysing the time evolution of atomic co- 
ordinates and velocities. Diffusion coefficients have been 
calculated using the asymptotic relation <| Ar(t) | 2 >= 
6Dt; i.e. by measuring the slope of the mean square 
displacement of the atoms as a function of time. Our 
results are displayed in table I for r s = 0.5,r s = 1 and 
r s = 1.2, in reduced plasma units D* = D/a 2 aj p i, where 

Upi{r s ) = (3/Af/) 1 / 2 r s 3 ^ 2 is the bare ionic plasma fre- 
quency. The results are shown in a log- log plot in Fig. 10. 
The relationship between D* and T = l/{r s T) follows 
quite accurately a power-law of the type D* = D T a . 
We have fitted our data to such an expression, obtain- 
ing for r s = 1 the following values: D a = (10.4 ± 1.4) 
and a = (—1.38 ± 0.07), and for r s = 0.5 the values: 
D = (4.0 ± 1.5) and a = (-1.37 ± 0.09). 

It is interesting to note that the diffusion coefficient 
follows the same relationship as in the OCP model. The 
OCP values for the parameters, fitted to classical MD 
simulations |5l| are D a = 2.95 and a = —1.33. The 
value of the exponent seems to be unaffected by elec- 
tronic screening, at least within the accuracy of our cal- 
culations. The prefactor, however, is clearly enhanced 
from its OCP value, and this can be readily understood 
in terms of the response of the electronic component to 
the motion of the protons. A rigid uniform electronic 
background (as in the OCP) does not have any influence 
on the dynamics of the protons. A polarizable back- 
ground weakens the proton-proton interaction thus in- 
creasing the mobility of the ions. In fact, the results 
obtained at r s — 0.5 are quite close to the OCP values, 
represented by the dashed line in Fig. 10. The differ- 
ence becomes much larger at r s = 1, where D* differs 



from its OCP values by a factor of 3; this contrasts with 
the Thomas- Fermi MD results of Zerah et al. f^fl , who 
observed a much milder effect (a factor 1.4 at T = 50 
and r s = 1). Therefore, the diffusion coefficient appears 
to be very sensitive to the treatment of electron screen- 
ing. The case of r s — 1.2 is slightly different, because 
the exponent appears to be larger than the OCP value. 
However, it can be seen that the error bars are also com- 
patible with an OCP-like power-law, represented by a 
line parallel to that of the other densities. 

The diffusion coefficient is also related to the ion veloc- 
ity autocorrelation function (VACF), which, in the OCP 
limit, exhibits a striking oscillatory behavior due to a 
strong coupling of the single-particle motion to the col- 
lective ionic plasma oscillations fill . Such oscillations 
were recently shown to persist at finite r s , by MD simu- 
lations using the approximate Thomas-Fermi kinetic en- 
ergy functional instead of the Kohn-Sham version fl52"| . 
The present ab initio calculations qualitatively confirm 
this behavior. In Fig. 11 we present the VACF for a 
typical simulation in the solid phase (T = 300 K) and 
then for three different temperatures in the fluid phase. 
The latter were computed in the supercell containing 
162 atoms. Finite size effects are not very significant 
as regards the general features of the VACF. However, a 
small frequency shift is observed, and decorrelation hap- 
pens faster in the larger sample. It is interesting to note 
that the fastest oscillation, i.e. the one associated with 
the ion plasma oscillations, is essentially temperature- 
independent. 

The power spectra of the ionic VACF are plotted in 
Fig. 12 for r s = 1, at several temperatures. The spectra 
exhibit a high-frequency peak (or shoulder at the highest 
temperature) at a frequency which amounts to 55 % of 
the bare ion plasma frequency (uj p i). The power spec- 
tra for r s = 0.5 are similar, with the difference that the 
high-frequency peak occurs now at a value which is 70 % 
of the bare plasma frequency p4| . As expected, electron 
screening shifts the vibrational spectra to lower frequen- 
cies. Temperature, however, does not affect the position 
of the high-frequency peak, which remains the only well- 
defined feature at high temperatures, while the rest of the 
spectrum merges into a structureless continuum. These 
spectra were obtained from 54 H atom simulations. The 
162 H atom sample yields a frequency 5 % lower than the 
54 H atom one. 



VIII. COLLECTIVE MODES: A SIGNATURE OF 
THE METAL-INSULATOR TRANSITION 

The oscillations in the VACF point to a long-lived lon- 
gitudinal collective mode, related to the ionic plasmon 
mode of the OCP pi] . We have computed the charge 
autocorrelation function: 

F zz (k,t) = -I < Pz(M) M-M) > (8) 
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where the Fourier components of the microscopic charge 
density are: 



Pz{k,t) = Pi(k,t) - p e (Kt) 



(9) 



with pi and p e the AIMD-generated time-dependent den- 
sities. In keeping with the Born-Oppenheimer approxi- 
mation, p e (k.,t) is a Fourier component of the expecta- 
tion value of the electronic density for the instantaneous 
ion configuration. In practice, we computed the average 
of Fzz (k, t) over the shell of equal- modulus k- vectors 



F ZZ (k,t) = J2 FzzM 



(10) 



From this, we computed the dynamical structure fac- 
tor Szz{k,u)) by Fourier transforming Fzz{k,t). In the 
r s = (OCP) limit, where the electrons form a uni- 
form (non polarizable) background, Szz(k^)i which 
there reduces to the dynamical structure factor of the 
bare ions, is simply the k-dependent spectrum of the 
ionic plasma oscillations |51j] . In the long wavelength 
limit the mode is undamped, and its characteristic fre- 
quency is the ion plasma frequency lo p \. Adiabatic elec- 
tron polarization transforms this ionic plasmon (or op- 
tic) mode into an acoustic mode for any finite value of 
r s This mode is to be identified with the famil- 

iar low-frequency ion-acoustic mode. Only if the system 
were treated as a fully dynamical ion-electron plasma, 
would the high-frequency plasma oscillation mode ap- 
pear, which is related to the fast electronic motions. 
This mode is obviously not accessible by adiabatic MD 
simulations. The conjectured scenario is confirmed by 
the results of our AIMD simulations. The dynamical 
structure factor Szz{k,w) was computed for the smallest 
wavenumber k compatible with the PBC (ka = 1.031 for 
the 54 atom system) and for selected larger wavenum- 
bers (ka < 3). The resulting Szz(k,oj) for r s = 0.5, 
r s = 1 and r s — 1.2 are shown in Fig. 13, for the 
smallest available wavenumber, and at a temperature 
just above melting. The sharp peaks are characteris- 
tic of the long-lived (weakly damped) mode anticipated 
above. The peaks shift to lower frequencies as r s in- 
creases due to enhanced electron screening, and in ac- 
cord with the behavior of the plasmon-like peak in the 
spectrum of the VACF (the single-particle excitation cou- 
pled to the collective plasmon mode) . The k-dependence 
of the spectrum Szz(k,ui) is illustrated in Fig. 14, for 
r s = 1; the resulting dispersion curve is shown in fig. 
15. A striking feature is the nearly constant width of 
the resonance peaks for ka < 1.5, pointing to a nearly 
k-independent damping mechanism. The damping in- 
creases dramatically at larger wavenumbers, while the 
dispersion curve bends over; the behaviour is reminiscent 
of that observed for a classical fluid of atoms interacting 
via an effective Yukawa potential |54|] . The bending over 
may be regarded as a remnant of the negative dispersion 
of the plasmon mode observed in the strongly- coupled 



OCP pi] . From the initial slope of the dispersion re- 
lation, we estimate a sound velocity of c s ~ 70 km/s 
at r s = 1, which is is consistent with the extrapolation 
of very recent results by Alavi et al. to the ultra-high 
1 
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Sound velocities are relevant to the 
obal free oscillations of Jovian planets, 
which have been recently measured for Jupiter [ p5[ . 

The remarkable feature is that this collective mode is 
much sharper than the usual sound mode observed in 
metals at comparable wavenumbers. Moreover, the peaks 
do not shift significantly with temperature, although they 
broaden. However, particularly for low values of k, the 
signature of the collective mode can be detectable up to 
quite high temperatures (well above 3000 K). This is 
seen in Fig. 16, where we plot the dynamical structure 
factor for ka = 1.238 (i.e. just before the bending-down 
point) in the 162-atom sample at r s = 1 and for three dif- 
ferent temperatures. It is interesting to note that this is 
precisely the temperature range where a transition is ex- 
pected to occur to the molecular (H 2 ) fluid phase at lower 
densities pSfl . The observed collective behavior may be 
regarded as characteristic of the metallic phase of hy- 
drogen, and is expected to change dramatically at the 
transition towards the molecular phase, which begins to 
show up at low temperatures at r s ps 1.3. Thus, an anal- 
ysis of Szz(k,uj) may provide an efficient diagnostic to 
locate the plasma phase transition at finite temperature. 



IX. CONCLUSIONS 

The main conclusions to be drawn from the present 
AIMD simulations of the hydrogen plasma in the high- 
density (r s < 1.2) regime may be summarized as follows: 

a) Due to the significant spacing between the quantized 
electronic states in the vicinity of the Fermi surface, the 
N-dependence of the statistical averages must be treated 
with great care, in order to extract meaningful results. 

b) A linear-response treatment of the ion-electron cor- 
relations yields reasonable results at r s = 0.5, but be- 
comes rapidly unreliable at lower densities. 

c) The bec structure, which is the stable low tempera- 
ture solid phase at least up to r s = 0.5, becomes unsta- 
ble at lower densities, where hep and simple-hexagonal 
phases appear. More work is needed to determine the 
full low-temperature phase diagram, also including zero- 
point-motion effects. 

d) The melting temperature drops sharply with de- 
creasing density, due to the enhanced efficiency of elec- 
tron screening of the effective interaction between ions. 
New interesting physics is likely to arise in the region of 
r s w 1.1 and T ss 100 — 200 K, where the existence of 
a bcc-hcp(fcc)-liquid triple point is argued, in a region 
where quantum effects in the protons are non-negligible. 

e) The fluid metallic phase behaves very much like a 
simple atomic liquid from a structural point of view, but 
the longitudinal collective dynamics of the ions retain a 
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strong plasma-like character at intermediate wavenum- 
bers. This reflects itself in unusually sharp peaks in the 
charge-fluctuation spectrum, which are gradually shifted 
to lower frequencies with decreasing density, as a result 
of electron screening. The damping, however, appears 
to be surprisingly insensitive to density, but is signifi- 
cantly enhanced by temperature. The acoustic charac- 
ter of the longitudinal mode is recovered at sufficiently 
small wavenumbers (feet < 1), in qualitative agreement 
with a simple linear screening picture. A strong damp- 
ing of the mode at intermediate wavenumbers should be 
a clear-cut signature of the plasma-to-molecular phase 
transition, which is expected to start at T — around 
r s = 1.3, and to move to finite temperatures of the order 
of a few thousand K at lower densities (r s > 1.3) 

f ) The single-particle motion of the ions couples to the 
longitudinal collective mode, and reflects itself in a strik- 
ing oscillatory behaviour of the velocity autocorrelation 
function, which is reminiscent of the behavior of the OCP, 
despite the action of strong electron screening. The re- 
sulting ionic self-diffusion constant is strongly enhanced 
at lower densities, for identical values of the plasma cou- 
pling constant T, but follows a power law similar to that 
observed in the OCP. 

The present AIMD simulations will be extended to 
lower densities, in order to characterize the plasma-to- 
molecular phase transition, starting from the high den- 
sity, metallic side. 
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FIG. 1. Electron distribution around the proton site for 
H in the perfect bcc structure at r 3 = 0.5, as a function of 
the plane wave energy cutoff. Dotted, dash-dotted and solid 
curves correspond to cutoffs of 60, 230 and 420 Ry, respec- 
tively. 

FIG. 2. Electron distribution around the proton site for H 
in the perfect bcc structure at r s = 0.5 (lower curves) and 
r s = 1 (upper curves), according to LDA (solid lines) and 
within linear response (dashed lines). Distances expressed in 
units of the ion-sphere radius a — r s as- 

FIG. 3. Electron distribution around the proton site for H 
at T=2000 K according to LDA (solid line) and LRT (dashed 
line), for r s = 1 (lower curves) and r s = 1.2 (upper curves). 



FIG. 4. Effective pair potentials for the proton-proton in- 
teraction at r s = 1: «lda fitted to AIMD trajectories (solid 
line) , and «lrt obtained within the RPA with local field cor- 
rections (dashed line). 

FIG. 5. Proton-proton pair distribution function at 
T — 1000 K and r 3 = 1 for different system sizes: (a) 54 
atoms (solid line), 128 atoms (dashed line), and 162 atoms 
(dotted line); (b) 128 atoms considering explicitly the whole 
open shell (81 states), where the lowest 64 are initially oc- 
cupied (solid line), 128 atoms considering explicitly only 
the 64 lowest occupied states (dashed line), and 162 atoms 
(dot-dashed line). 

FIG. 6. Fermi-Dirac distribution function for T = 1000 K 
(solid line), T = 5000 K (long-dashed line), and T = 10000 K 
(dot-dashed line). The electronic density of states associated 
with the last two occupied shells is also shown (short-dashed 
line), as well as the single shot of the Kohn-Sham eigenval- 
ues during an AIMD simulation at T = 1000 K (vertical 
bars). The EDOS is broadened with a Gaussian window 1 eV 
wide. 

FIG. 7. Proton-electron pair correlation function at fi- 
nite temperature (solid and fluid phases) at r 3 — 0.5 and 
T = 1000,2500,5000 and lOOOOif (lower curves), and at 
r s = 1 and T = 200,500,1900 and 73007^ (upper curves). 
The lower curve of each set corresponds to the lower temper- 
ature. 

FIG. 8. (a) Proton-proton pair distribution function at 
r s = 1 and T = 300, 500, 1000, 1500, 2000, 2500 and 7300 
K. (b) Integral of the above distribution function for the 
same temperatures. The function N(r) indicates how many 
protons are contained in a sphere of radius r, on average. 

FIG. 9. Proton-proton pair distribution function at r 3 = 1 
for a 162 H atom sample, as a function of temperature: 
T = 1000 K (solid line), T = 2000 K (dashed line) and 
T = 3000 K (dot-dashed line). 



FIG. 10. Reduced diffusion coefficient at r s — 0.5 (circles), 
r s = 1 (triangles), and r s = 1.2 (diamonds) in a log-log plot. 
Solid lines are the results of a linear regression fit on the 
logarithms. The dashed line is the OCP result. 

FIG. 11. Velocity autocorrelation function at r s = 1 for 
T = 300 K (dotted line), T = 1000 K (dashed line), 
T = 2000 K (dot-dashed line), and T = 3000 K (solid line). 
The abscisa (time) is expressed in units of the plasma fre- 
quency w p i. 
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FIG. 12. Power spectra of the velocity autocorrelation 
function at r s — 1, obtained by FT the velocity autocorre- 
lation function. The curve at T = 300 K (solid line) corre- 
sponds to the solid phase, while the remaining curves cor- 
respond to T = 800 K (dotted), 1000 K (short-dashed), 
2000 K (long-dashed), 2500 K (dot-dashed) and 7300 K 
(solid). These are all in the fluid phase. Frequencies are in 
units of the ionic plasma frequency u) p \. Curves are normal- 
ized to integrate to the number of ionic degrees of freedom. 

FIG. 13. Charge-charge dynamical structure factor at the 
smallest available wave vector, ka — 1.031, for three values of 
r s - 

FIG. 14. Charge-charge dynamical structure factor at 
r s = 1 for different values of momentum: ka = 0.715 
(solid line), 1.011 (dotted), 1.238 (short-dashed), 1.599 
(long-dashed), 2.022 (dot-dashed), 2.476 (diamonds), and 
2.948 (crosses), from a simulation with 162 H atoms. The 
curves have been obtained by Fourier transforming the auto- 
correlation functions (8), filtering with an appropriate decay- 
ing exponential those corresponding to the largest wavevec- 
tors in order to reduce noise. All curves are normalized by 
the static structure factor Szzik) = J Q Szz(k,cu) dui, such 
that the area under every curve is 1. 

FIG. 15. Dispersion relation for the collective ion acoustic 
mode in Szz{k,u)), at r s — 1. Filled circles are points corre- 
sponding to the k- vectors allowed by the PBC, in a simulation 
box containing 162 H atoms. The solid line is a guide to the 
eye, arising from a 4 th -degree polynomial fit to the data. 

FIG. 16. Dynamical structure factor at r s — 1 for 
ka = 1.238 at T = 1000 K, (dot-dashed line), T = 2000 K 
(solid line), and T = 3000 K (dashed line). 



TABLE I. Diffusion coefficient as a function of the plasma 
coupling parameter F. The error in the determination of the 
diffusion coefficients is AD* = 0.001. 
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